Sampling Chaotic Trajectories Quickly in Parallel 

J. MachtaQ 

Department of Physics, University of Massachusetts, Amherst, MA 01003-3720 

(Dated: February 8, 2008) 

Abstract 

The parallel computational complexity of the quadratic map is studied. A parallel algorithm is 
described that generates typical pseudotrajectories of length t in a time that scales as logt and 
increases slowly in the accuracy demanded of the pseudotrajectory. Long pseudotrajectories are 
created in parallel by putting together many short pseudotrajectories; Monte Carlo procedures 
are used to eliminate the discontinuities between these short pseudotrajectories and then suitably 
randomize the resulting long pseudotrajectory. Numerical simulations are presented that show the 
scaling properties of the parallel algorithm. The existence of the fast parallel algorithm provides a 
way to formalize the intuitive notion that chaotic systems do not generate complex histories. 
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I. INTRODUCTION 



Chaotic systems cannot be predicted for very long times because of the exponential 
divergence of nearby trajectories. Associated with the divergence of trajectories is a lack 
of history dependence; the current behavior of the system is not dependent on the past 
behavior. The absence of history dependence can be understood in various ways. Here I 
take a computational perspective on chaotic systems and analyze trajectories in terms of 
computational complexity. If, using a massively parallel computer, a typical long trajectory 
can be manufactured in far fewer parallel steps than the actual length of the trajectory then 
the trajectory lacks history dependence or historical complexity. Alternatively, if parallelism 
does not allow one to generate a typical trajectory much more quickly than its actual length 
then the trajectory displays a complex history dependence. 

These considerations are illustrated using the one-dimensional quadratic map, 

x n+l = rx n (l - x n ) = f(x n ) (1) 

with x n G [0, 1] and < r < 4. The objective is to produce typical trajectories for the 
map but since computing devices are necessarily restricted to finite precision we are really 
interested in generating typical pseudotrajectories JTJ. A pseudotrajectory of accuracy 5 is 
a sequence {y n \n — 0, • • • , t} such that for all n (0 < n < t), 

\y n+1 - f(y n )\<S. (2) 

We suppose that the parallel computer works to a precision substantially better than 5 so 
that bounds such as given in Eq. |2| can be checked with reasonable certainty. 

For chaotic dynamics with a positive Lyapunov exponent A, a pseudotrajectory and an 
exact trajectory that are initially equal remain close only for a time that is roughly given by 
(log<5)/A. Nonetheless, for small 5 a pseudotrajectory will have nearly the same statistical 
properties as a real trajectory and, in any case, numerical results about chaotic systems are 
learned from pseudotrajectories not real trajectories. 

The goal then is to produce a pseudotrajectory chosen from the uniform distribution 
over pseudotrajectories. A typical pseudotrajectory of length t can be generated using one 
processor in time linear in t by iterating the map using arithmetic of precision much better 
than 5 and then adding a noise term on each step chosen from the uniform distribution on 
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[—5,5]. Could a typical trajectory be produced in far fewer than t parallel steps with the 
aid of many processor? 

The model of parallel computation implicit in this question is the PRAM (parallel random 
access machine), the standard model in the theory of parallel computational complexity [F2J. 
A PRAM is an idealized, fully scalable device with many identical (except for distinct 
integer labels) processors. Processor all run the same program and all communicate with 
a global memory in unit time. Massive parallelism is envisioned here however the number 
of processors is required to be polynomially bounded in tlog(l/5), the effective number of 
(binary) degrees of freedom of a pseudotrajectory of length t and accuracy 5. 

In the next section we show how to produce a typical trajectory in parallel using Monte 
Carlo path sampling. The procedure is correct but inefficient. We then describe how sim- 
ulated annealing together with path sampling can produce a typical pseudotrajectory in 
parallel time that scales linearly in logt and polynomially in log(l/<5) using a number of 
processors that is polynomial in tlog(l/5). 



II. PATH SAMPLING OF PSEUDOTRAJECTORIES 

The method described in this section is based on path sampling ideas put forward by 
Chandler and collaborators ||. The uniform probability density for pseudotrajectories , 
PG/0,2/1, ■■■,Vt) is given by 

t 

V(y , 1/1, ... , y t ) = p(yo) II p (yn\y n -i) (3) 

n=l 

where 

P Wy) J-^-^< S (4, 
I otherwise 

and p(y) is the invariant distribution. 

A simple Monte Carlo procedure can be used to sample paths from V. Given an existing 
pseudotrajectory, a single time n > is chosen and a proposal for a new value for y n is 
obtained according to 

y'n = f{Vn-i) + e (5) 

where e is chosen from the uniform distribution on [—5, 5} . The proposed value is accepted 
if it is also the case that \y n+ i — f(y' n )\ < 8. It is straightfoward to verify that this Monte 
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Carlo procedure satisfies detailed balance with respect to V. The question of ergodicity of 
the Markov chain in the space pseudotrajectories is less clear. However, even if ergodicity 
holds, the actual mixing time for the Monte Carlo procedure would be long when 8 is small 
since the time to obtain an independent pseudotrajectory is at least as great as 1/S 2 , the 
time to diffuse a distance order one given a step size of 5. To be considered an efficient 
process for generating pseudotrajectories , the parallel time should increase no more rapidly 
than some power of the logarithm of 1/5. This goal can be achieved by first using a simulated 
annealing procedure to produce a pseudotrajectory that is random on long time scales then 
to use the above Monte Carlo path sampling to randomize the small scales. 

III. SIMULATED ANNEALING FOR PSEUDOTRAJECTORIES 

Long pseudotrajectories are constructed by independently generating many segments or 
short trajectories and then "welding" the segments together. In the welding step the dis- 
continuity between successive segments is eliminated by simulated annealing. In addition to 
simulated annealing, it is sometimes necessary to extend segments to obtain a weld to the 
next segment. 

The fundamental time scale in the system is 

r = -(log*)/|A|, (6) 

the time required for typical errors to grow to be order one. The length of the segments, K 
used in the construction should be longer than the fundamental time r so that the beginning 
and end of each segment is uncorrelated. We also need to allow for extensions at either end 
of the segment and for a "warm-up" so that the initial point of the segment is chosen from 
invariant distribution. Thus, in practice, for each segment we start with a random number 
and iterate the map L > K times, choosing the segment of length K from a predetermined 
part of the longer sequence of length L. 

Having made a collection of segments, we now attempt to weld them together into a long 
pseudotrajectory . This is done in such a way that the initial value or anchor point of each 
segment is held fixed. The discontinuity between successive segments is annealed until all 
errors are less than 5. The Monte Carlo annealing procedure is designed to lower the error 
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e(y n -i,y n ,y n+ i) associated with three successive elements, y n -i,y n and y n +i, 

e(yn-i,y n ,yn+i) = \y n - f{yn-i)\ + \y n +i - f{y n )\- (7) 

If e(y n -i,y n ,y n+ i) < 5 for every n then we have a pseudotrajectory. For each time n, with 
the exception of the anchor points, the Monte Carlo annealing procedure begins with a 
measurement of e = e(y n -x, y n , y n +i)- If e < 5 nothing is done. Otherwise a new value y' n is 
proposed, 

y'n =Vn + e (8) 

where e is chosen as a Gaussian random variable with mean zero and standard deviation 
e/2. If e' = e(y n -i,y' n ,y n+ i), is less than e, the proposal is accepted as the new value for 
y n . If the error increases, the proposal is accepted with probability e _ ^ e '~ e ) . The value of 
the inverse temperature for each Monte Carlo step and is taken to be j3 = e/2 so that the 
acceptance ratio is independent of the size of the error. In a single Monte Carlo sweep, all 
except the initial and final values of each segment are processed using the above procedure. 
Given t/2 processors this can be done in constant parallel time by first processing the even 
and then the odd values of the time n. 

In the chaotic regime (A > 0) the annealing procedure should yield a valid pseudotrajec- 
tory for large enough K and sufficiently many Monte Carlo sweeps. In practice, however, 
some welds require a very large number of sweeps. Specifically, the probability distribution 
for the number of sweeps needed to achieve a weld has a long tail leading to parallel running 
times for creating pseudotrajectories that are dominated by the few most difficult welds. 
Two additional kinds of steps, called forward shifts and backward shifts cure this difficulty. 
Suppose segment s together with its final condition, the first element of segment s + 1, is 
not fully annealed after a predetermined number of annealing sweeps. Then segment s is 
restored to its original state and either it is extended forward one step or segment s + 1 is 
extended backwards one step. In the case of a backwards shift, the element prepended to 
segment s + 1 is considered a new anchor points and serves as the new final condition for 
segment s. The net effect of either a forward or backward shift is that the discontinuity 
between segments s and s + 1 occurs with a different pair of numbers. Annealing sweeps and 
shifts are interleaved, a fixed number of Monte Carlo annealing sweeps are attempted and if 
all errors are not less than 5, a shift is done. The process is repeated until a satisfactory weld 
is achieved. Successive shifts are alternately of the forward and backward type. In Sec. [V], I 
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show that the combination of Monte Carlo annealing and shifts produces a pseudotrajectory 
in O(logi) parallel steps. 

Shifts serve several purpose. First, they simply provide for the possibility of more Monte 
Carlo sweeps, though if this were their only function it could be accomplished by directly 
increasing the number of sweeps. Second, shifts permit the algorithm to perform properly 
for periodic orbits or nearly periodic stretches of aperiodic orbits. The annealing procedure 
by itself cannot generate long pseudotrajectories for periodic orbits since welds are often 
attempted between segments that are out of phase with one another. Adding shifts to the 
annealing procedure insures that periodic pseudotrajectories will be correctly generated. For 
example, consider the case of a period two orbit, a single forward or backward shift of some 
segments will insure that all welds are satisfactory. For period d orbits, as many as d — 1 
shifts are necessary to insure that all welds are satisfactory. 

Shifts may also provide padding around hard to weld regions of a trajectory. During a 
shift, new points are added to one end of a segment but no points are removed. Thus shifts 
do not bias the pseudotrajectory against difficult to weld regions in the invariant measure. 
For example, it is observed that if the final condition for a segment is very near the maximum 
of the support in the invariant measure at r/4 then one or more backward shifts are usually 
necessary so that the point near r/4 is surrounded by a region of small errors and is not 
involved in the annealing process. 

IV. FULL PARALLEL ALGORITHM FOR PSEUDOTRAJECTORIES 

This section provides the details of the parallel algorithm for producing typical pseudo- 
trajectories that combines the path sampling method of Sec. y and the annealing procedure 
of Sec. |ITT[ First the annealing procedure is used to generate a pseudotrajectory and then 
path sampling is used to further randomize it. The algorithm is controlled by several param- 
eters: t is the total length of the desired pseudotrajectory , 5 is the desired accuracy, K is 
the length of each segment, K' is the warm-up length, E, an even number, is the maximum 
number of shifts that are attempted, M\ is the number of Monte Carlo annealing sweeps 
carried out between shifts and Mi is the number of Monte Carlo path sampling sweeps. The 
algorithm is described below: 

1. In parallel, generate S — \t/K~\ sequences {a;W} each of length L = K' + K + E, 

6 



with s = 0, • • • , S — 1 and m = 0, • • • , L — 1. The initial value of each sequence is a 
uniform random number on (0, 1) and subsequent values are obtained by iterating the 
map L — 1 times to precision much greater than 5. This step requires O(L) parallel 
time. 

2. These S sequences are used to define ES segments {yn' q ^} each of length K, where 
the index q, < q < E gives the number of shifts. y£' q ) = Xn+E/2+K'+\q/2\ f° r 
< n < K — 1 while y^'-i — xS e/2+k> -\q/2]+K> ■ Note that the final point in segment 
y(s,g) j s a taken from sequence s + 1. 

3. In parallel, for each s < S, and each q < E, anneal segment y^ s ' q \ The annealing 
procedure consists of Mi Monte Carlo sweeps. During a single annealing sweep first 
the even and then the odd elements of the segment are updated in parallel. The initial 
and final points y^'^ and y^'-i are held fixed during the annealing procedure. A single 
Monte Carlo update of the point y^'^ consists of the following procedure: 

• Compute e = e(y^l\,y^ q \y^\) from Eq. 0. 

• If e < S, do nothing. If e > 5 propose a new value y' = y^ q "> + e where e is a 
Gaussian random variable with mean zero and standard deviation e/2. 

• Compute e' = e(y^_ 1 , y', ym+i)- If e' < e accept the proposed change, y^ q ^ <— y'. 
lie' > e accept the proposed move with probability exp[— /3(e'— e)] where (3 — e/2. 

The parallel time required for this step is O(Mi). 

4. In parallel, for each s < S find Q(s), that smallest value of q such that the segment 
is successfully annealed. The annealing is successful for this segment if, for all < 
m < K — 1, the errors are sufficiently small, e{y^^ y^ q \ < ^- If f° r an y s i 
annealing is unsuccessful for all q < E, the algorithm fails. This step can be carried 
out in constant parallel time. 

5. The full pseudotrajectory y* is a concatenation of sequences obtained from each orig- 
inal sequence x^ s \ The contribution to the pseudotrajectory from is the con- 
catenation of {x$\E/2 + K' - \Q(s - l)/2] < m < E/2 + K' + [Q(s)/2\} and 
{Vm® |0 < m < K — 1}. The first of these sequences is composed of the anchor 
points and the second sequence is the annealed segment. To obtain a pseudotrajectory 
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of length exactly t, the pseudotrajectory obtained above is simply truncated after t 
steps. 

6. The path sampling Monte Carlo procedure described in Sec. [TT] further randomizes 
y* m . During a single sweep, first all the even and then all the odd elements of the 
pseudotrajectory are updated in parallel. The number of sweeps is M 2 . On each 
Monte Carlo step a new value for y* m is proposed according to Eq. || and accepted only 
if the trajectory is still a pseudotrajectory within error 5. The randomization step 
requires parallel time 0(M 2 ). 

V. VALIDITY AND COMPLEXITY OF THE PARALLEL ALGORITHM 

The central questions addressed in this section are (1) whether the algorithm succeeds in 
creating a pseudotrajectory (2) how the scaling of the number of parallel steps depends on 
the length and accuracy of the pseudotrajectory and the parameter r of the map and (3) 
whether the algorithm samples the uniform distribution on pseudotrajectories. A sequential 
algorithm that carries out the annealing and path sampling routines one segment at a time 
was used to study these questions. In the simulations reported below, the parameters are 
chosen to be Mi = M 2 = 5r 2 , K = 5r and K' = 1000. The assumption behind these choices 
is that memory is lost on a time scale r so that placing independently chosen anchor points 
separated by K = 5r is satisfactory. The annealing process that welds successive segments 
is expected to influence a region whose length is order r. Since information is transmitted 
diffusively by local Monte Carlo moves, having the number of Monte Carlo sweeps scale as 
t 2 should suffice. 

The annealing stage of the parallel algorithm can be studied one segment at a time since 
each segment is independently annealed. First, I observed that, given enough shifts, the 
annealing step always produced a successful weld. The choice of the maximum number of 
shifts E for the annealing stage must be large enough to make the failure probability for 
the whole algorithm small. For long trajectories, t ^> r, the choice of E is determined by 
the tail of the distribution of the number of shifts, Q required to obtain a weld since the 
whole procedure fails if even one segment is not successfully annealed. Suppose C(-) is the 
cumulative probability distribution for Q. An estimate of the maximum number of shifts, 
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FIG. 1: Logarithm of the complement of the cumulative distribution for the number of shifts, 
log 10 (l — C(Q)) vs. number of shifts, Q. From left to right, 5 = 10 -5 to 1CT 10 . 

E needed to insure all segments are annealed is given by the relation 

(t/K)(l-C(E))<l. (9) 

Figure]]] shows log 10 (l — C(Q)) vs. Q. These data were collected for the case of r = 3.7 (the 
period doubling transition to chaos occurs at r = 3.5699 . . .) and the six curves from left to 
right are for 5 = 10~ 5 through 10~ 10 , respectively. Each curve is obtained from annealing 
10 5 segments except for the 5 = 10~ 10 curve which is obtained from 6 x 10 4 segments. For 
r = 3.7 the Lyapunov exponent is A = 0.354 and so, for example, with 5 = 10~ 7 , r = 45.5, 
K = 228 and Mi = 10366. Over a reasonable range following an initial transient and before 
the noise becomes large, the data falls on straight lines suggesting that the distribution is 
asymptotically exponential, C(Q) ~ 1— exp(— Q/cr). Equation ^ then implies that E ~ a Int 
and we can conclude that the parallel running time is O(logt) since no other contribution 
to the running time depends on the overall length of the pseudotrajectory. 

How does the decay constant a and thus the running time depend on the choice of the 
accuracy 5. Figure fj shows o vs. 5 on a logarithmic scale for the case r = 3.7 and suggests 
that a is a polynomial function of log 5. The other simulation parameters, L, Mi and M 2 
are also polynomial in log 5 so we conclude that the full algorithm has a running time that 
is polynomial in log 5 and linear in logt. 

I also considered two other parameter values for the quadratic map, r = 3.6 where A = 
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FIG. 2: Decay constant a for the distribution of shifts vs. accuracy 5 on a logarithmic scale. 

0.183 and r = 3.95 where A = 0.577. In both cases, the accuracy was set to 5 = 10~ 7 . The 
decay of (1 — C(Q)) appears to be exponential in both cases but with rather different values 
of the decay constant: o = 450,15 and 9.5 for r = 3.6,3.7 and 3.95, respectively. Either 
a depends strongly on A or perhaps there are additional r dependent factors controlling a. 
For example, for r = 3.7 and 3.8 the invariant measure has support on a single interval but 
for r = 3.6 the support consists of two intervals. 

The annealing stage of the algorithm creates a pseudotrajectory but it is not typical in 
the sense of being chosen from the distribution of Eq. |3|. On long time scales, the pseudotra- 
jectory is randomized by the random choice of initial conditions for each sequence. However, 
the individual errors, e n = — f(Un) ar e not guaranteed to be independent random vari- 
ables on the interval [—5, 5}. For example, anchor regions of the pseudotrajectory have errors 
much less than 5. The hypothesis is that M 2 = 0(r 2 ) path sampling Monte Carlo sweeps are 
sufficient to randomize the short time scales and produce a typical pseudotrajectory from 
the pseudotrajectory produced by the annealing stage. To check this hypothesis, I computed 
mean values and autocorrelation functions for errors and cross correlations between errors 
and values of the pseudotrajectory . The quantities (e n ), ((e 2 ) —5 2 /3), (e n+ ie n ), (e n+ i?/*) 
and ((^n+i e n) ~ ^ 4 /9) were all found to be zero within error bars for the case r = 3.7 and 
5 = 10 -7 . Here the angled brackets indicate an average over segments and over n. The van- 
ishing of these quantities is a necessary but not sufficient condition that pseudotrajectory 
is chosen from the uniform distribution described by Eqs. |3] and |j. More work is needed to 
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firmly establish that the algorithm with 0(t 2 ) path sampling sweeps samples the uniform 
distribution to good approximation. However, even without the path sampling stage, the 
pseudotrajectories produced by the annealing stage are typical in a different sense. As shown 
in PL any pseudotrajectory shadows an exact trajectory (i. e. remains close to over its entire 
length) though possibly for a larger value of r. 

VI. CONCLUSIONS 

I have exhibited a parallel algorithm that generates pseudotrajectories of the quadratic 
map. Numerical evidence suggests that the parallel time required to generate a typical 
pseudotrajectory increases linearly in logt and polynomially in \og(l/S) though more work 
would be required to establish these scalings with certainty. Essentially the same parallel 
algorithm can be applied to other one- dimensional and higher dimensional maps. It would 
be interesting to explore whether the annealing/shift procedure is sufficient to efficiently 
sample pseudotrajectories for other maps. 

Since there is little demand for very long pseudotrajectories of the quadratic map, the fast 
parallel algorithm is probably not of practical value. The significance of its existence and 
complexity is that it characterizes the history dependence of the map. The existence of a fast 
parallel simulation is a strong statement against history dependence since it shows that the 
logical path from independent random numbers (used to drive the Monte Carlo procedures) 
to a typical pseudotrajectory is much shorter than the length of the pseudotrajectory. The 
length of this logical path is one measure of the potential for generating historical complexity. 
In very few logical steps, very little complexity can be arise. The idea that complexity tends 
to emerge at the "edge of chaos" j|] is born out here since the basic time scale r for the 
parallel algorithm diverges when the Lyapunov exponent vanishes. An appealing feature 
of characterizing systems by computational complexity is that very different systems in 
statistical physics systems, for example diffusion limited aggregation |J, sandpiles |J or the 
Bak-Sneppen model |7J|, can be compared to one another within the same framework. 
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